% This script estimates parameters for each generator in a loop

format longE

global forward searchhistory

% Estimation scripts and calls
%  estimateNew  (data handling and looping over generators)
%       -MC.m  = calculates MC given the heatrate and marginal costs
%       -estimate_policy.m = starts estimation routine (sets starting values for the value function Vlast)
%           -binify2 = binifies prices into non uniform bins possibly changing the number of bins
%            -EstTransMat3,2, or 1.m = different algorithms for estimate the price transistion matrix; EP specifies which one
%            -fminsearchbig(electML_Policy) = search to maximize the likelihood
%                 -electML_Policy
%                     -electllfi_policy  %change code here to have decision be implemented now or 1 period in the future (foward)
%                         -V_solve_sf  % matrix inversion is done here
%                             -Profit_Policy
%                        
%             -llvar(electllfi_policy) = calculate standard errors


% Speed: i={1,2}
%   laptop: i=1: 127 secs; i=2: 261 secs
%   desktop nomultithread: i=1: 83 secs i=2: 134secs

%%  SET SPECIFICATION
% choose the firm number i=1,..,n, if not looping
i=22;   %Witchita
i=14;  %PB2SES
i = 1; % Calenergy

% forward algorithm where this period's decision is not encated until next period (1-yes, 0=no)
forward = 0;
policy = 1;
spec = 'linear_norm_fixed'


% compare solution times for policy vs value function in solutiontimes.xls

% hourly or not
hourly = 1;
months = 3;
EP = 1;   %use EstTransmat_linear (1), EstTransMat2 (2), EstTransMat3 (3),EstTransMat4 (4), EstTransMat_nonparam(0)

% Size parameters
P = 40;  %NOTE: this is not the final number of price bins. Binify2.m may reduce the number of price bins
if hourly
  disp 'Hourly'
  I = 24;
  beta = 0.999994;    %for .95 annual with hourly data
  %beta=0.5
  % beta = 0.1
  % I can get the value functions to converge to e-5 to e-7 with this
   % discount factor
else
  disp 'Interval'
  I=96;
  beta = 0.9999985;    % for .95 annual with 15 min data
  % I can get the value functions to converge to e-4 to e-5 with this
   % discount factor (with hourly data)
end
S = 2;


% can get policy function iteration down to e-11 with discount factor .99
% can get policy function iteration down to e-12 with discount factor .9



%file = 'C:\Documents and Settings\Joseph\My Documents\Econ\Environ\Data\ERCOT\';
file ='';
% Import parameter estimates for all units for the given month
%       Parameters for each unit should be in the following form
%           Name, MC, LSL, HSL, START, SIGMA
%[Name, MC,LSL,HSL,START,SIGMA] = textread('unitParameters.csv', '%s %f %f %f %f %f ');

%% READ IN DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % Operating Parameters: TX subscripted parameters come from a private
 % company. Data shared by Steve Puller of TX A&M
[name,LSL,HSL,LOL,HOL,RampRate,StartTime]=textread('./Data/Operate.csv', '%s %f %f %f %f %f %f ');
 %[name,LSL,HSL,LOL,HOL,RampRate,StartTime,LSL_TX,HSL_TX,RampRate_TX]=textread('./Data/Operate2.csv', '%s %f %f %f %f %f %f %f %f %f ');

% read fuel prices
[month_fuel,year_fuel,gas_cost,coal_cost]=textread('./Data/Fuel.csv', '%f %f %f %f ');
f=[month_fuel,year_fuel,gas_cost,coal_cost];
% read heat and emissions characteristics
%[name,HR_EPA,HR_EIA,NoxR ,So2R,Co2R_EPA,Co2R_EIA,gas]=textread('./Data/Rate.csv', '%s %f %f %f %f %f %f %f')

[name,HR_EPA,HR_EIA,NoxR ,So2R,Co2R_EPA,Co2R_EIA,gas, HR_TX, So2R_TX, NoxR_TX , Co2R_TX ]=textread('./Data/Rate2.csv', '%s %f %f %f %f %f %f %f %f %f %f %f');

%calculate MC for each generator
% use TX data
HR = HR_EIA; %HR_TX
co2ratio = Co2R_EIA./HR_EIA;
Co2R = co2ratio.*HR; 
temp = [Co2R, Co2R_EIA, Co2R_EPA];
VOM=0;

% use parameters created by get_params.do which takes the "best parameters
% from  CEMS, Steve's data, and EIA.
%[name,HR2,Co2R2, LSL2, HSL2, VOM, gas]=textread('./Data/newparams.csv', '%s %f %f %f %f %f %f ')
%[name,HR,Co2R, LSL, HSL, VOM, gas]=textread('./Data/newparams.csv', '%s %f %f %f %f %f %f ')




% name and estimating month 
[name,month_est]=textread('./Data/month_est.csv', '%s %f ');
N = length(name)

%% SET UP FOR ESTIMATION  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% month of data used to estimate costs for each generator
notestimated = (month_est==0);
keepmonth = month_est ;  
keepmonth(notestimated)=1;
if months==3
    disp '3 month'
    fuelmonth=5;    %use fuel prices for the beginning month. over the three month period
elseif months==1
    disp 'single month'
    fuelmonth=keepmonth;   %a vector of monts
end


% set parameters for emissions costs
So2cost=600;  %$/ton
Noxcost= 2000; %$/ton
Co2cost =0;   %$/ton

% create vector of fuel prices
FuelCost = gas.*gas_cost(fuelmonth)+(1-gas).*coal_cost(fuelmonth);



MC = calcMC(HR, NoxR, So2R, Co2R, Noxcost, So2cost, Co2cost, FuelCost, VOM);
% MC2 = calcMC(HR2, NoxR, So2R, Co2R2, Noxcost, So2cost, Co2cost, FuelCost, VOM);
% 
% [MC,MC2,LSL,LSL2,HSL,HSL2]


%number of generators
gens = size(MC,1);

% a vector of starting values
sStart= ones(gens,1)*160680;
sSigma = ones(gens,1)*21000;
sFC = ones(gens,1)*-6160;
%initial_params_v=[sStart,sSigma, sFC];
initial_params_v=[sStart,sSigma, sFC];
%initial_params_v=[sStart,sSigma, sFC, MC];

% create file name as characters as opposed to cells

if hourly
   % for HOURLY
    
    strname = strcat('./Data/ercotdata_', name, '_h.csv'); 
else
    % for INTERVAL
    strname = strcat('./Data/ercotdata_', name, '.csv');
end
strname=char(strname);

%% Estimate for firm i 

% a hold for results
num_params = size(initial_params_v,2);
allresults=zeros(num_params,2,N);


%matlabpool local8 8  %for UA cluster
%parfor i = 1:N
%for i = 2:N    
%for i = [1,12,15,21,22]
for i =  1:N
%for i=1
    i
    
    if i==2
        continue
    end
    
    if month_est(i)==0
        % if I am not estimating this one than skip it
       continue
    end    
        
    %current file
    inputfilename = strtrim(strname(i,:))
    
    %'./Data/'
    
    %open dataset (interval)
    [price, interval, state, action, month, year, timeid] = textread(inputfilename, '%f %f %f %f %f %f %f');    
    %[price, interval, state, action, month, year, timeid] = textread('ercotdata_WITCHITA.csv', '%f %f %f %f %f %f %f');    
  
  
    % What is this all about?? 9/15/2010 (can't figure it out so I am commentingout)
    % change action and state to reflect forward decision process
  
    
    
    % which months to use for estimation
    if months==1
        % for a single month
        keep = (month==month_est(i)) .* (year==2006);
    elseif months==3
        % for three months
        keep = (month==5 | month==6 | month==7 ) .* (year==2006);
    else
        quit
    end
    
    if forward==1
        disp 'forward algorithm deployed'  %see logbook 12/15/2011 for explanation
        state = action(1:end-1);  %state is now the observed operating condition for the unit which was decided by last periods action
        action = action(2:end); %the action is now the observed operating state next period  which is when the action is implemented c
        price = price(1:end-1);  % the rest of the variables need to be shortened so the vectors are the same size
        interval = interval(1:end-1);
        month = month(1:end-1);
        year = year(1:end-1);
        timeid = timeid(1:end-1);
    end
    
    
    % select the subset of data for the months selected above 
    indicies = find(keep); %get indicies of  non-zero elements
    price = price(indicies);
    interval = interval(indicies);
    state = state(indicies);
    action = action(indicies);
    timeid = timeid(indicies);
    month = month(indicies);
    year = year(indicies);
   
    X = [price, interval, state, action, month, year, timeid];
 
    %estimate parameters
    initial_params = initial_params_v(i,:)';
   
    
    % output to AMPL file
    %outputAMPL(name(i),price, interval, state, action, beta, MC(i), LSL(i), HSL(i), P,I,S, initial_params);
 
    % display stats for current estimator
    genname = name(i)
    mc=MC(i)
  
    if hourly
        % for HOURLY
        lsl=LSL(i)
        hsl=HSL(i)
        %lsl=LSL_TX(i)
        %hsl=HSL_TX(i)
    else
        % for INTERVAL
        lsl=LSL(i)/4;
        hsl=HSL(i)/4;
    end

    tic
    results = estimate_policy(strname(i,:),price, interval, state, action, beta,mc,lsl, hsl, P,I,S, initial_params, EP, forward, policy)
    toc 
    allresults(:,:,i)=results;
    
    % save the results (saves all the matrices in the work space)
    if months==1 
        file = strcat('est_', name(i,:),'_',num2str(i),'_1month_',date,spec,'.mat')
        filename=char(file)
        save(filename, 'results')
    elseif months==3
        file = strcat('est_', name(i,:),'_',num2str(i),'_3months_',date,spec,'.mat')
        filename=char(file)
        save(filename, 'results')
    end
    
end
 
   if months==1 
        file = strcat('allresults_1months',date,spec, '.mat')
        save(file)
    elseif months==3
        file = strcat('allresults_3months',date,spec, '.mat')
        save(file)
        
    end



